%% Clean workspace
    clc
    clear all

%% Load parameters
% sample_manager
    manager
    freqlisttemp=[];
    for i=1:hn
        freqlisttemp=[freqlisttemp, freqlist*i];
    end
    freqlist=freqlisttemp;
    perlist = 1./freqlist;
%%
    datelist = startdate:timestep:enddate;

    % resulting spectrogram
    if ifgpu
        res=gpuArray(ones(1,length(perlist),9));
        resc=gpuArray(ones(1,length(perlist),2,9));
    else
        res=(ones(1,length(perlist),9));
        resc=(ones(1,length(perlist),2,9));
    end

    stdlist=ones(length(datelist),9);
    meanlist=ones(length(datelist),9);
    lN=ones(length(datelist),9);
    baddays=ones(length(datelist),9)*NaN;
    baddays(:,1)=mjuliandate(datevec(datelist'));
    maini=1;
           
%% Event list
    % add new fields to the end of the structure!!!
    event=struct('objID',{},...
                 'MJD',{},...
                 'freq',{},...
                 'freqer',{},...
                 'okno',{},...
                 'range',{},...
                 'detector',{},...
                 'profile',{},...
                 'rms',{});
    eventnum=0;
%% main cycle

maini=maini-1;
tic

while maini<length(datelist)
    maini=maini+1;
    try
    %% reading data
    t=datelist(maini);
    bool=0;
    for i=floor(t):1:ceil(t+anint)
        if bool==0
            mas=readbinkw(i, datadir);
            bool=bool+1;
        else
            temp=readbinkw(i, datadir);
            temp(:,1)=temp(:,1)+86400000*bool;
            mas=[mas; temp];
            bool=bool+1;
        end
    end
    
    % cut beginning and end
    st   = (t-floor(t))*24*3600*1000;
    temp = abs(mas(:,1)-st);
    stp  = find(temp==min(temp));
    if stp<1
        stp=1;
    end
    
    et   = (t+anint-floor(t))*24*3600*1000;
    temp = abs(mas(:,1)-et);
    etp  = find(temp==min(temp));
    if etp>length(mas)
        etp=length(mas);
    end
    
    mas=mas(stp:etp,:);
	
	% Move timeframe
    mas(:,1)=mas(:,1)-mean(mas(:,1));
	
	% Convert into seconds (originally in ms)
	mas(:,1)=mas(:,1)/1000;
	
	% Convert into flux (counts/s)
    for j=2:9
		mas(2:end,j)=mas(2:end,j)./(mas(2:end,1)-mas(1:end-1,1));
	end
	mas=mas(2:end,:);
	
	% Close all files (important!)
    fclose('all');
    
    %% Mark good data
    good=mas*0;
    meanmas=mean(mas);
    trendmas=mas;
    for i=okna
        [ trendmas(:,i), good(:,i), baddays(maini,okna) ] = TrendDef(mas(:,1),mas(:,i));
    end
    
    %%
    for okno=okna
        %%
        goodness=(good(:,okno)==1);
        data= mas(goodness,okno);
        data=data-trendmas(goodness,okno);
        stdlist(maini,okno)=std(data);
        meanlist(maini,okno)=mean(trendmas(goodness,okno));
        lN(maini,okno)=length(data);
        %% do analysis?
        if doanalysis
            if ifgpu
                times=gpuArray(mas(goodness,1));
                data=gpuArray(data./stdlist(maini,okno));
            else
                times=(mas(goodness,1));
                data=(data./stdlist(maini,okno));
            end
            peri=5000;
            for i=1:peri:length(perlist)
                tempi=i:min(length(perlist),i+peri-1);
                [tempa, tempb]=FourierMod(times,...
                                       data,...
                                       perlist(tempi),...
                                       ifgpu);
                res(1,tempi,okno)=1/2*(tempa+tempb);
    %             resc(1,tempi,1,okno)=tempa;
    %             resc(1,tempi,2,okno)=tempb;
            end
            is=SourceMarker(freqlist,gather(res(1,:,okno)),lN(maini,okno),0);
            temp=struct2cell(is);
            tempid=temp(2,:,:);
            tempid=cell2mat(tempid(:));
            for id=1:100
                indexes=find(tempid==id);
                if ~isempty(indexes)
                    %%
                    er=[];
                    freq=[];
                    for ind=indexes(:)'
                         harmnum=is(ind).hn;
                         magnit=is(ind).magn;
                         er=[er pi./sqrt(3)./(anint*24*3600)/harmnum/sqrt(magnit)];
                         freq=[freq is(ind).freq/harmnum];
                    end
                    eventnum=eventnum+1;
                    event(eventnum).objID=id;
                    event(eventnum).MJD=mjuliandate(datevec(datelist(maini)))+anint/2;    
                    event(eventnum).freq=sum(freq./(er.^2))/sum(1./(er.^2));
                    event(eventnum).freqer=(sum(er.^-2))^(-0.5); 
                    event(eventnum).okno=okno;   
                    event(eventnum).detector=(okno>5)+1;
                    event(eventnum).range=okno-1-(event(eventnum).detector-1)*4;   
                    %% Profile extraction
                    event(eventnum).profile=EP(times,data.*stdlist(maini,okno),event(eventnum).freq,32);
                    event(eventnum).rms=std(event(eventnum).profile(:,1));
                end
            end
        end
    end    
    %%
    toc
    disp(maini/length(datelist))
    end
end

